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The Kohn-Sham orbital kinetic energy density r a (r) = JT Wi a | V^i ff (r) | is one fundamental 
quantity for constructing meta-generalized gradient approximations (meta-GGA) for use by density 
functional theory. We present a computational scheme of r CT (r) for full-potential linearized aug¬ 
mented plane wave method. Our scheme is highly accurate and efficient and easy to implement to 
existing computer code. To illustrate its performance, we construct the Becke-Johnson meta-GGA 
exchange potentials for Be, Ne, Mg, Ar, Ca, Zn, Kr, Cd atoms which are in very good agreement 
with the original results. For bulk solids, we construct the Tran-Blaha modified Becke-Johnson po¬ 
tential (mBJ) and confirm its capability to calculate band gaps, with the reported bad convergence 
of the mBJ potential being substantially improved. The present computational scheme of r CT (r) 
should also be valuable for developing other meta-GGA’s in FLAPW as well as in similar methods 
utilizing atom centered basis functions. 

PACS numbers: 71.15.Mb, 77.22.Ej 


I. INTRODUCTION 


Nowadays, density functional theory has become the 
dominant method for electronic structure calculations. 
Accompanying its 50 year development, the quality of 
the approximated energy functionals are also constantly 
improved. 2 In order of higher accuracy, Perdew has cat¬ 
egorized various functionals in terms of the “Jacob’s 
ladder” P The first rung of the ladder is the local den¬ 
sity approximation (LDA) which depends only on the 
spin densities (p|, p^). The second rung consists of gener¬ 
alized gradient approximations (GGA) which depend on 
(pf, pi, Vp-f, VpjJ. The third rung consists of “meta gen¬ 
eralized gradient approximations (meta-GGA)” which 
further include the Kohn-Sham orbital kinetic energy 
density T a (r), so that the functionals now depend on 
(Pt,Ph V A-Vp4.,r t ,T4.). 

As a fundamental component of meta-GGA, r CT (r) 
should be implemented in a way which is highly accu¬ 
rate and efficient at the same time. For the full-potential 
linearized augmented planewave method-'"' 6 a previous 
implementation of r CT (r) existf 7 . By definition r a (r) re¬ 
quires the gradients of all occupied Kohn-Sham orbital 
(w is the occupation number): 


Tcr( r) = ^ w i0 

2=1 


V^ i(7 (r) 


(1) 


Since direct computation of the gradients on dense 
meshes in the unit cell would be too costly, the authors 
make use of the Kohn-Sham equation 

V 2 1 

+v C r(r)>ip ia (r) = e ia ip ia (r) (2) 

to convert T CT (r) to the following form 

Arfi) = ^V 2 p CT (r) - v a {r)p a (r) - X e ia p ia (v) (3) 

i—1 


which involves only the total spin density p (7 (r), the 
Kohn-Sham effective potential u CT (r) and the Kohn-Sham 
orbital eigenvalues e,^. Eq.([3]) is the form being imple¬ 
mented. It is more efficient than Eq.([T]) since the depen¬ 
dence on individual orbital has been removed. 

Nevertheless, it is often unnoticed that the conver¬ 
sion from 0 to 0 is not always valid, since typically 
not all electrons states are solved by Kohn-Sham equa¬ 
tion. In fact, core states are usually solved with rel¬ 
ativistic corrections which is the case in the standard 
FLAPW implementation. Then, expression 0 is even 
not non-negative-definite. Invalid numeric values are 
mostly found around nuclei where the Dirac-Kohn-Sham 
solution deviates the most from the Kohn-Sham solution. 
Technically, these invalid data may be zero’ed without af¬ 
fecting the band structure very much, for which the dom¬ 
inant effects come from potential in the chemical bond¬ 
ing regions. It is unavoidable, however, that the conver¬ 
gence of self-consistency is often ruined, and some author 
questions whether the bad convergence is intrinsic to the 
meta-GGA exchange potential itselfP 

Of course, one can intentionally use Kohn-Sham equa¬ 
tion to solve not only the valence states but also the core 
states. With the sacrifice of the relativistic corrections, 
Eq.0 becomes equivalent to Eq.(Jl]). Even so, however, 
the equivalence is purely mathematical rather than nu¬ 
merical. This is because in almost all cases Kohn-Sham 
equation is solved only approximately by use of finite 
basis sets. Through the conversion from 0 to ([3| the 
incomplete basis set error is brought into the values of 
T er( r )- 

It is always desirable to calculate r cr (r) directly from 
the original definition Eq. 0 since this equation is valid 
regardless of the details of how the orbital are solved. 
The purpose of this work is to present a computational 
scheme of r a {r) for full-potential linearized augmented 
planewave method (FLAPW). Our scheme is highly ac- 
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curate and efficient and requires minimum change to the 
existing computer code. Moreover, the increased accu¬ 
racy of r a (v) substantially improves the convergence of 
self-consistency. 

This paper is organized as the following: In section II 
the FLAPW method is briefly reviewed. As this method 
has become very popular in recent years and there have 
been numerous literatures about itp our introduction is 
limited to the minimum content which is relevant to the 
present work. In section III, we derive the formulae of 
our computational scheme of r CT (r) for both the valence 
and the core states. The performance of this scheme is 
illustrated in section IV. We first construct the Becke- 
Johnson meta-GGA exchange potential for atoms, and 
then the Tran-Blaha modified version of the potential 3 
for bulk materials, by which we discuss its major feature 
for calculating semiconductor band gaps. We summarize 


our work in Section V. Within the FLAPW framework, 
the unit cell is partitioned into atomic spheres (called the 
muffin-tin region) and the space in between (called the 
interstitial region). Implementation of T a (r) in the inter¬ 
stitial region is simple and straightforward. Therefore, 
most discussions are focused on the muffin-tin region. 
Throughout this work atomic units are used. To sim¬ 
plify the notations, spin index and some other symbols 
are suppressed wherever confusion is avoided. 

II. BRIEF INTRODUCTION TO THE FLAPW 
METHOD 

The FLAPW method solves the valence states and 
core states by different techniques. Valence states are 
expanded by the following basis functions: 


4( r ) 
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lm 
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ai(g)ui(r) +bi(g)ui(r ) 
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(muffin-tins) 

(interstitial) 


(4) 


Each basis function is of a hybrid form: Inside the muffin- 
tins it takes the linear combination of atomic orbital, 
while within the interstitial it is simply plane wave. In 
this fashion, the basis function best accounts for the gen¬ 
eral feature of the wave functions in the whole space 
which behave like atomic wave functions close to nuclei 
and become free-electron-like in-between. The optimal 
shape of the basis function allows the expansion of the 
valence state 

^nk(r) = 5> nk (G)0 g (r), (g = k + G) (5) 

G 

to use small cutoff for the reciprocal lattice vector G. 

In Q , Y lm is spherical harmonics which is widely used 
in atom centered basis functions. The ui(r) is solved by 
the radial differential equation with a pre-chosen energy 
parameter and keeping only the spherical component of 
the potential. The combination it;(r)l) m (f) is called a 
muffin-tin function. The other radial function ui(r) is 
the energy derivative of ui (r). It serves as the first order 
correction to u;(r), so that the basis function becomes 
applicable to a wide energy range. The use of the two 
muffin-tin functions inside one atomic sphere is proposed 
by Anderson 5 as the linearised approximation to the orig¬ 
inal augmented plane wave method of Slatei®. The linear 
coefficients ai( g) and bi( g) are so chosen that the ba¬ 
sis function and its first derivative are both continuous 
across the sphere boundary. 

To account for space group symmetry, density and po¬ 
tential are expanded by symmetrized plane waves (called 
star functions) in the interstitial region and symmetrized 
spherical harmonics (called lattice harmonics) in the 


muffin-tin region. In particular, within an atomic sphere, 
a lattice harmonics is constructed as 

K v (r) = '^2,c v {m v )Y lvm „{i:) ( 6 ) 

m u 

so that it is invariant under any local point group sym¬ 
metry IZi: 

HiK v (r) = K v { f) (7) 

Expand the density and the potential by the lattice har¬ 
monics 

P( r ) = '52P-'( r ) K A*) ( 8 ) 

v{r) = ^v v {r)K v {r) (9) 

symmetry requirements are automatically fulfilled. 

By applying IZi to Eq. (|3j) , we realize that r(r) has the 
same symmetry as p{ r) and v(r), and therefore can also 
be expanded by the same set of lattice harmonics 

T ( r ) = ^2 T A r ) K A^) ( 10 ) 

V 

The expansion coefficients are calculated by: 

T„(r) = J r(r)K*(r)dn 

Core states in the standard FLAPW implementation 
are solved by Dirac-Kohn-Sham equation: 

{ca-p + c 2 /3 + w (0) (r)}^ KjU (r) = ( n ) 
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to better account for relativistic effects. In 0 , a and 
/3 are the standard Dirac matrices, v ^ (r) only contains 
the spherical component of u(r). The four component 
solution to Eq.(ll) is 




9i K {r)x K ^{i) 

ifi K {r)v r XKiJ,(*) 


( 12 ) 


III. IMPLEMENTATION SCHEME OF r(r) FOR 
THE FLAPW METHOD 

A. r”(r) of the muffin-tin region: contribution from 
the valence states 

Within the muffin-tin region, the expansion of the va¬ 
lence state can be written as: 


with g nK (r) and / nre (r) being the radial part of the major 
and minor components, respectively. Note that the core 
states are assumed to be completely restricted within the 
muffin-tins so that their wave functions have no intersti¬ 
tial part. Besides, <? nK (r) and f nK {r) are not expanded 
by basis set, but are directly integrated on the real space 
grid. 

The total kinetic energy density is contributed by the 
core states and all occupied valence states: 

r(r) = r c (r) + T v (r) (13) 


fpn k(r) = Y 


Ai m (nk)ui{r) + Bi m {nk)ui(r) 


Yi m { r) 


in which we have defined the coefficients 


Ai m (nk) = Y ^nk(g) a/(g) Y^(g) 

G 

B lm (nk) = Y z nk(s) Mg) Y * m ( g) 

G 

The total kinetic energy density of the valence states in¬ 
side the muffin-tin region can be written as: 


= Y w ” k v 0nk(r) 


= VV( r)~Y w n 

nk 


V£k( r )VVnk(r) - ^ nk (r)V 2 V’* k (r) 


(15) 


Note that to reach at ( flbj ) we did not make use of any On the right hand side of Eq.(15), the first term is 
physics equation such as Kohn-Sham equation. There- already calculated by GGA. The remaining two terms 

fore, 0 is equivalent to (|T|) regardless of the details of can be derived by use of the following relation: 
how the orbital are solved. 


l [w{r)Y lm { r)} = (w"(r) + 


2 W'{r) l(l + l)W(r) 


Yim{ r) 


(16) 


Eq.( 161 is a favorable property of all muffin-tin functions, 


by which the laplacian operation only affects the radial 
part while leaves the angular part unaltered. Using (16) 


for p~5| ), we get t v ( r) in the muffin-tin region. Projected 
to the lattice harmonics representation, the final expres¬ 
sion of the expansion coefficients is: 
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For comparison, we also write out the expansion coefficients of the density: 


'v! £ 

'-T 

II 

UiUl> 

^ J ^nk ^ ^ A-Vm' -A-lm i^lm 
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UlUV 

nk mm' 

^ ^ ^nk ^ ^ ^l'm' ^lm i^lm 

Yl v m„ 

+ 

UlUv 

nk mm' 

^ ^ ^nk ^ ^ m' ^lm i^lm 

Y lvmu 

+ 

UlUv 

nk mm' 
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(18) 


Eqs.(fl7|) and (18) are similar. Their differences are lim¬ 
ited to the radial parts while their angular parts are the 
same. This is not a coincidence, but a consequence of 
the property (16) of the muffin-tin functions. Eq. © 


is 


the expression being implemented in the present work. 
It is highly accurate because it is equivalent to Eq. 0 
and no extra error is introduced throughout the deriva¬ 
tion. Moreover, the expressoin of the lattice harmonics 
expansion also ensures the correct symmetry of r(r). 

Eq.© is equally efficient because r^(r) can be cal¬ 
culated together with pj)(r). During the construction of 
the valence density, almost all computational cost of (18) 


is spent on the angular part, while the radial part takes 
no more than a few percent. By calculating r"(r) to¬ 
gether with p„(r), the costly angular part needs to be 
evaluated only once. The extra cost for the radial part 
of ( p~7] ) is very small. Note that the radial part of ( fl7| ) 
has no orbital dependence since it only requires the ra¬ 
dial basis functions ui{r) and iii{r) which are universal 
for all valence states. Technically, the implementation of 
the radial part of ( [T7| to the construction of the valence 
density is also straightforward and requires little human 


labor. 


B. T c (r) of the muffin-tin region: contribution from 
the core states 

Core states also contribute to the r(r) of the muffin-tin 
region. Although core states are in the four-component 
form, each row is a separate muffin-tin function so that 


Eq.(16) is equally applicable. At ground state, all core 
sub-shells are fully filled. Consequently, core density is 
spherical and r(r) also is because they must have the 
same angular parts: 


r c (r) = r 0 c (r) 


(19) 


The expansion coefficient of the first lattice harmonics 
A'o(r) is contributed by all occupied sub-shells: 


)(0 =© r o,i( r ) 


( 20 ) 


For the l = k sub-shell: 


Ur) = 


21 

\/~U 


, 2 l{l + 1 ) 2 , 2 l{l — 1 ) f2 


9iK + fL + 


( 21 ) 


For the l = —n — 1 sub-shell: 


c / \ 2 ( l + 1 ) J / 2 1(1 + 1) 2 if' 2 , + 1 )(^ + 2 ) f2 

r 0 ,i(r) — s T ^ 9in Y Jm T ^2 /m 


( 22 ) 


Since the radial solution of Dirac-Kohn-Sham equation 
is slightly irregular, r(r) diverges at the nucleus which 
is similar to that of all GCA potentials. Alternatively, 
core states can also be solved by Kohn-Sham equation. 
Then the expression of Tq (r) can be directly derived from 


the valence expression © in which one only needs to set 
= 0 and l = l 1 and remove the linearization term iij (r). 

C. r v ( r) of the interstitial region 

In the interstitial region, only valence states contribute 
to r(r): 
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r v (r) = ^E w ^E( k + G ‘) ( k + G^z^iGMGj) e^~ G *> r 

L nk GiGj 


(23) 


Eqs.@, @ and ( |23[ ) have been implemented 
to the FLAPW code of Northwestern University.^ To test 
the performance of our computational scheme, we next 
construct the Becke-Johnson type meta-GGA potentials 
for atoms and bulk semiconductors. For atoms, core 
states are solved by Kohn-Sham equation for compar¬ 
ison with the earlier work. For semiconductors, both 
Kohn-Sham equation and Dirac-Kohn-Sham equation 
are tested for solving the core states, and we find that 
the results as well as the speed of convergence are nearly 
the same. Therefore, throughout this work all core states 
are solved by Kohn-Sham equation. 


IV. THE BECKE-JOHNSON TYPE META-GGA 
POTENTIALS 


A. The atomic potentials 


The Becke-Johnson (BJ06) potential is a meta-GGA 
exchange potential attempting to approach the optimized 
effective potential (OEP) 12 of atoms. BJ06 consists of 
two terms: 


<f 6 (r)=<r(r) 



(24) 


The first term is the Becke-Roussel exchange potential 13 
which is derived by assuming hydrogen-like exchange 
hole. u® 306 (r) is determined by a nonlinear equation 
involving p CT (r), Vp CT (r), V 2 p CT (r) and r a ( r), and closely 
resembles the Slater averaged exchange.^ Compared to 
OEP, u®^ 89 (r) is too deep and lacks the characteristic 
shell-structures. Therefore, the second term is added 
for correction. This term is repulsive which reduces the 
strength of the exchange. Besides, it strongly varies 
across atomic shells which restores the OEP shell struc¬ 
tures. Summing up the two terms, the BJ06 potential is 
close to the accuracy of the OEP of atoms, yet is much 
simpler since it is directly constructed from local quan¬ 
tities rather than being solved by the complicated OEP 
equation. Nevertheless, BJ06 is a model exchange poten¬ 
tial without corresponding exchange energy. Therefore, 
it cannot be used for calculating total energy. 

For Be, Ne, Mg, Ar, Ca, Zn, Kr, Cd atoms, the BJ06 
potentials have been published by Becke and Johnson. 10 
We have generated the same potentials by our FLAPW 
code and the results are plotted in Figure [l] As can be 
seen, both the location and the height of the “bumps” be¬ 
tween atomic shells are well reproduced. Besides, close to 
the nuclei all potentials have the correct, asymptotically 
flat shape, indicating the adequacy of the treatment of 
r(r) in the muffin-tin region. The depth of the potential 




FIG. 1. Becke-Johnson exchange potential for Be, Ne, Mg, 
Ar, Ca, Zn, Kr, Cd atoms calculated by the FLAPW method 
with the implementation of t„ (r) of this work. 


bottom also agrees well with the original results. Over¬ 
all, all atomic potential are in excellent agreement with 
RefP. The remaining slight differences may be ascribed 
to the use of different boundary conditions in the two 
groups of calculations: While the original calculations are 
performed with the molecular code, the present calcula¬ 
tions use the FLAPW code with supercells and periodic 
boundary condition. 


B. The bulk potentials 


To test our implementation of r(r) for bulk materi¬ 
als, we construct the Tran-Blaha modified version of the 
BJ06 potential (mBJ) 11 for solids. The TB09 exchange 
is essentially the same as BJ06, with only the weights of 
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the BR89 term and the shell term being changed: 


TB09 

"x,cr 


(r) = 


^ R89 (r) + (3c - 2) 1 J A (25) 

7T V 12 V p a (r) 


In Eq. (25) the weights are determined by the c parame¬ 


ter, which is essentially adjustable, but can also be cal¬ 
culated “pseudo- ab initio- ly” by empirical formulae pro¬ 
posed by the authors®. For c = 1 TB09 reverts to the 
original BJ06. The total mBJ potential is formed either 
by TB09 exchange alone, or by combining it with a corre¬ 
lation potential such as LDA. It is such slight adjustment 
of the weights which leads to the surprising discovery that 
the mBJ potential is capable of calculating semiconduc¬ 
tor band gaps. Albeit being a local potential so that its 
computational cost is essentially maintained at the LDA 
level, in many cases the mBJ potential can achieve band 
gap accuracy which is even comparable to the much more 
sophisticated GW approximation. 

Applications of the mBJ potential has been quickly 
growing in recent years which also motivates the present 
work. However, the existing FLAPW implementation 
uses Eq.(3]) to calculate r(r), and therefore suffers from 
the problems mentioned in the Introduction. It is thus 
desirable to recheck the mBJ potentials by our imple¬ 
mentation of r(r). In all our calculations, the c parame¬ 
ter of J25| is automatically determined by the suggested 
formula™ 


c = A + Bg 


(26) 


with A = 0.488, B = 0.5, and 


- = J_ f 1 < W(r)| 
9 Kell J 2 ^ pt( r ) 

cell 


]V/A(r)|\ 

pH r) J 


d 3 r 


(27) 


For the correlation potential we use LDA. 16 One differ¬ 
ence with the earlier work is that in our calculations spin- 
orbit coupling (SOC) is included perturbatively through 
the second variational approach. Typically, the effect of 
SOC on the value of band gap is smaller than 0.1 eV. But 
for materials (such as ZnTe) with heavy elements SOC 
can change the bandgap by as much as 0.3 eV. 

In Figures [2] and [3] we plot 1 ^ the spin-up kinetic en¬ 
ergy density and the mBJ potential of the NiO (001) 
plane for comparison with the similar plot of Ref® At 
the nuclei r(r) shows sharp spikes while u mBJ (r) achieves 
the minimum, although neither quantities diverges. The 
present implementation of r(r) allows us to reveal even 
the finest details of the mBJ potential in Figure [3] Espe¬ 
cially, all shell structures are clearly seen. Another fea¬ 
ture is that, although r(r) and u mBJ (r) are very smooth 
within both the muffin-tin and the interstitial regions, 
across the sphere boundary they are obviously discontin¬ 
uous. Such discontinuity is also observed in the potential 
plot of NiO in Ref.® and is irrelevant to the potential 
itself or to the implementation of r(r). Rather, it is 
caused by the discontinuity of the basis functions. By 
Q , the FLAPW basis function is continuous only to the 
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O- Ni 


*1 • 
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FIG. 2. (Color online) Rubber-sheet plot of the up-spin ki¬ 
netic energy density log 10 Tj- of the NiO (001) plane. Atomic 
positions are illustrated at the left bottom of the plot. 





FIG. 3. (Color online) Rubber-sheet plot of the up-spin mBJ 
potential u) nBJ of the NiO (001) plane. Atomic positions are 
illustrated at the left bottom of the plot. 


first derivative while the meta-GGA require the second 
derivative of the density. In fact, the basis function is 
already discontinuous even at the zero’th order because 
in Eq. @ the summation is cutoff at finite (/to). By the 
same reason, the GGA potential is also discontinuous 
across sphere boundary, although we notice that usually 
the discontinuity of the meta-GGA is more serious than 
the GGA. 

With our highly accurate r(r), we have calculated the 
mBJ band gaps for all semiconductors in the SC40 test 
set.® The results are plotted in Figure [J] together with 
the LDA band gaps to compare with experiment. Indeed, 
for all semiconductors mBJ systematically improves the 
band gaps. Especially, the improvement is nearly as good 
for the smallest band gap (of InSb with E g = 0.23 eV) as 
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FIG. 4. (Color online) Comparison between the mBJ band 
gaps, the LDA band gaps and the experimental band gaps for 
the 40 semiconductors of the SC40 test set. 

for the largest (of MgO with E g = 7.31 eV). Besides, we 
also confirm that in most cases the size for the band gap 
increases nronotonically with increasing c. Nevertheless, 
serious errors still exist for the three Ba compounds, BaS, 
BaSe, BaTe, suggesting that there might be problems in 
the Ba potential or the experimental data. 

We have thus found that the improvement of r(r) does 
not qualitatively change the most important feature of 
the mBJ potential. This is not surprising since the im¬ 
provement of r(r) is most significant around the nuclei 
while the band gap correction mainly depends on po¬ 
tential far awayf 8 In more details, the dominant contri¬ 
bution to the band gap correction comes from the shell 
term which is repulsive so that it pushes all band states 
upward. Away from the nuclei the shell term increases 
and for open systems it approaches a positive constant at 
r —> oo. Because the conduction states are usually more 
delocalized than the valence states, the upshift of the con¬ 
duction states by the shell term is larger than the valence 
states, and therefore the band gap is enlarged. Through 
this mechanism, band gap correction is not substantially 
affected by the improvement of r(r) around the nuclei. 
For the same reason, it is neither affected very much by 
the divergence of r(r) if core states are solved by Dirac- 
Kohn-Sham equation. 

It has been reported before that self-consistency by the 
mBJ potential is harder to achieve than LDA or GGAP 
This is confirmed in our work. Typically, using the mBJ 


potential 50% more iterations are needed to achieve self- 
consistency than LDA or GGA. For an example, for M 0 S 2 
LDA takes 22 iterations to converge, while mBJ with our 
implementation of r(r) requires 32 iterations. Neverthe¬ 
less, we have also tried implementing r(r) by Eq. (3]h by 
which we found that the speed of convergence is drasti¬ 
cally ruined. For M 0 S 2 the mBJ calculation with Eq. © 
requires 79 iterations to converge, while calculation of 
NiO does not even converge at all. In fact, in all our 
testings, the implementation of r(r) by ([l]) always con¬ 
verges better than by ©, and the speed of convergence 
is not sensitive to whether the core states are solved by 
Kohn-Sham equation or Dirac-Kohn-Sham equation. Be¬ 
sides, for more than 70 solids we have calculated with ([!]), 
all converge well except for CoO and FeO which fail be¬ 
cause mBJ is incapable of lifting the near degeneracy of 
the 3d states. Therefore, it is safe to conclude that at 
least a large part of the reported convergence problem is 
not intrinsic to the potential but is caused by technical 
reasons. 


V. SUMMARY 

In this work we have presented an highly accurate and 
efficient computational scheme for the Kohn-Sham or¬ 
bital kinetic energy density r CT (r) to the full-potential 
linearised augmented plane wave method. To test its 
performance, we have constructed the Becke-Johnson ex¬ 
change potential for atoms which are in very good agree¬ 
ment with the original results. For bulk solids we have 
constructed the Tran-Blaha modified Becke-Johnson po¬ 
tential for semiconductors and confirmed its capability 
to calculate semiconductor band gaps. As to the conver¬ 
gence problem reported before, we found that a large part 
is due to the impropriate implementation of r(r). With 
respect to accuracy, efficiency, easiness of implementa¬ 
tion and speed of convergence, our scheme all supersedes 
the existing implementation. We expect this work to be 
valuable for developing other meta-GGA’s in FLAPW as 
well as in similar methods utilizing atom centered basis 
functions. 
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